Загрузите датасет very_low_birthweight.RDS (лежит в папке домашнего задания). Это данные о 671 младенце с очень низкой массой тела (<1600 грамм), собранные в Duke University Medical Center доктором Майклом О’Ши c 1981 по 1987 г. Описание переменных см. здесь. Переменными исхода являются колонки ‘dead’, а также время от рождения до смерти или выписки (выводятся из ‘birth’ и ‘exit’. 7 пациентов были выписаны до рождения). Сделайте копию датасета, в которой удалите колонки с количеством пропусков больше 100, а затем удалите все строки с пропусками.
very_low_birthweight <- readRDS("very_low_birthweight.RDS")
full_cols <- very_low_birthweight %>% summarise_all(~ sum(is.na(.))) %>%
select (!c(lol, magsulf,meth,toc, pvh, ivh, ipe)) %>% names () # колонки, где меньше 100 пропусков
very_low_birthweight_с <- very_low_birthweight %>% select(all_of(full_cols)) %>%
mutate (id = seq(1,671)) %>% # введем колонку id для обозначения пациентов до всех удалений
filter (if_all (where(is.numeric),
~!is.na(.))) # удалим строки, где есть NA в численных значениях
# ранее использовала drop.na(), но так учитываюися и категориальные переменные (на 4 кейса меньше остается), а нам все же важнее complete cases по количественным переменным
#summary (very_low_birthweight_с)
#very_low_birthweight_с %>% tbl_summary(by = dead) %>% add_p()
very_low_birthweight_с %>%
skimr::skim()
| Name | Piped data |
| Number of rows | 536 |
| Number of columns | 20 |
| _______________________ | |
| Column type frequency: | |
| factor | 4 |
| numeric | 16 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| race | 3 | 0.99 | FALSE | 4 | bla: 304, whi: 212, nat: 13, ori: 4 |
| inout | 0 | 1.00 | FALSE | 2 | bor: 449, tra: 87 |
| delivery | 1 | 1.00 | FALSE | 2 | vag: 271, abd: 264 |
| sex | 1 | 1.00 | FALSE | 2 | mal: 268, fem: 267 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| birth | 0 | 1 | 84.63 | 1.54 | 81.51 | 83.43 | 84.77 | 85.82 | 87.48 | ▅▆▇▇▅ |
| exit | 0 | 1 | 84.76 | 1.55 | 81.05 | 83.56 | 84.89 | 85.99 | 87.72 | ▂▆▇▇▅ |
| hospstay | 0 | 1 | 46.76 | 63.28 | -295.00 | 21.00 | 40.00 | 64.00 | 797.00 | ▁▇▁▁▁ |
| lowph | 0 | 1 | 7.22 | 0.13 | 6.53 | 7.14 | 7.22 | 7.32 | 7.55 | ▁▁▃▇▂ |
| pltct | 0 | 1 | 204.07 | 80.70 | 16.00 | 146.00 | 204.00 | 256.00 | 571.00 | ▃▇▅▁▁ |
| bwt | 0 | 1 | 1137.02 | 239.76 | 400.00 | 967.50 | 1160.00 | 1331.25 | 1500.00 | ▁▃▆▇▇ |
| gest | 0 | 1 | 29.27 | 2.23 | 23.00 | 28.00 | 29.00 | 31.00 | 38.00 | ▂▇▇▁▁ |
| twn | 0 | 1 | 0.20 | 0.40 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▂ |
| apg1 | 0 | 1 | 5.01 | 2.65 | 0.00 | 2.00 | 5.00 | 7.00 | 9.00 | ▅▆▅▇▇ |
| vent | 0 | 1 | 0.54 | 0.50 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00 | ▇▁▁▁▇ |
| pneumo | 0 | 1 | 0.17 | 0.38 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▂ |
| pda | 0 | 1 | 0.20 | 0.40 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▂ |
| cld | 0 | 1 | 0.26 | 0.44 | 0.00 | 0.00 | 0.00 | 1.00 | 1.00 | ▇▁▁▁▃ |
| year | 0 | 1 | 84.63 | 1.54 | 81.51 | 83.44 | 84.77 | 85.82 | 87.48 | ▅▆▇▇▅ |
| dead | 0 | 1 | 0.12 | 0.33 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| id | 0 | 1 | 328.72 | 183.24 | 2.00 | 175.50 | 329.50 | 479.25 | 671.00 | ▇▇▇▇▆ |
#альтернативный способ расчитать дни в госпитале
#very_low_birthweight_с %>% mutate(term = as.numeric(date_decimal(exit)- date_decimal(birth))/86400)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
# t-test
stat.test <- compare_means(
lowph ~ inout, data = cleaned_data ,
method = "wilcox.test"
)
stat.test
## # A tibble: 1 × 8
## .y. group1 group2 p p.adj p.format p.signif method
## <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 lowph born at Duke transported 0.000000297 0.0000003 3e-07 **** Wilcox…
# Create a box plot with rstatix
ggplot(cleaned_data, aes(inout, lowph, color = inout)) + # This is the plot function
geom_boxplot(width = 0.2, fill = "white", alpha = 0.1) +
labs(x = "Method", fill = "Method") +
stat_pvalue_manual(
stat.test,
y.position = 7.7, step.increase = 1,
label = "p.adj"
)
Н0 - медианы (т.к тест Манн-Уитни) значений lowph для новорожденных,
рожденных в Дюке и перевезенных не различаюся. Н0 отвергаем ( p-value
очень мало).
Вариант интерпретации - дети транспортированные в госпиталь имеют более низкий показатель lowph (самый низкий pH крови в первые 4 дня жизни) и более низкую выживаемость, чем дети, рожденные в госпитале. Вероятно из других больниц чаще поступают дети в более тяжелом состоянии состоянии, чем рожденные в самом центре. Низкие значения pH крови характерны для тяжелых нарушений гомеостаза, что вероятно требует госпитализации в медицинский центр при университете.
require(psych)
require (corrplot)
corr_cl_data <- cleaned_data %>% select (where (is.numeric)) # новый датасет в количественными и ранговыми переменными
correlation <- psych::corr.test(corr_cl_data, method = "spearman", adjust = "fdr") # c поправкой на множ. сравнения
# график матрицы корреляций
corrplot(corr = correlation$r,
type="upper", order="hclust", method = "number", col = COL2('PiYG'),
p.mat = correlation$p, sig.level = 0.01, insig = "blank")
# network (показывает силу связей и взаимоположение)
library(corrr)
net <- corr_cl_data %>%
cor() %>%
network_plot()
net
btw и gest значительно скорелированы, hospstay и bwt/gest слабо
отрицательно скорелированы.
net plot визуализирует взамосвязи и степень корреляции переменных.
birth_scaled <- scale(corr_cl_data) # шкалирование
birth_dist <- dist(birth_scaled) # вычисление дистанций
df_dist.hc <- hclust(d = birth_dist,
method = "ward.D2")
fviz_dend(df_dist.hc,
cex = 0.6)
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
## Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
#library("ape")
#plot(as.phylo(df_dist.hc), type = "fan")
Не очень удачныый способ представления большого количества наблюдений.
library(pheatmap)
hierarchy <- pheatmap(birth_scaled,
show_rownames = FALSE,
clustering_distance_rows = birth_dist,
clustering_method = "ward.D2",
cutree_rows = 5, #зададим количество кластеров для рядов
cutree_cols = length(colnames(birth_scaled)),
angle_col = 45,
main = "Dendrograms for clustering rows and columns with heatmap")
hierarchy
На основе дистанций по средним шкалированных значений переменных
построили дендрограмму для переменных (сверху). Видно, что дистанция
между bwt и gest минимальна (по средним значениям), и эти переменные
“похожи” между собой и визуально представляют один кластер, далее скорее
всего в один кластер объединены переменные apg1, lowph и pltct.
Hostpstay отстоит от остальных переменных максимально далеко по
значениям матрицы дистанций.
По такому же принципу построили дендрограмму для строк, то есть иерархически объединили наблюдения (отдельных новорожденных) по степени близости на матрице дистанций. Значение кол-ва кластеров (отсечки) 5 задали самостоятельно в настройках.
PCA- метод снижения размерности, то есть преобразования данных таким образом, чтобы, снизив количество переменных, сохранить наибольшее число информации об “вариантивности” данных, то есть соблюсти отношения между строками/наблюдениями. Также метод помогает преобразовать данные таким образом, чтобы убрать корреляцию между переменными (в случае “birth” это найденные ранее скоррелированные переменные).
В результате PCA мы получаем “новые” нескоррелрованные переменные - компоненты и их количество равно кол-ву колонок. Компоненты имеют иерархию по признаку того, какой процент дисперсии исходных данных они объясняют (Proportion of variance). Метод позволяет анализировать и строки, и колонки единовременно.
Шкалирование нужно, чтобы корректно соотносить и делать математичекие операции разных переменных друг с другом.
library(factoextra)
library(FactoMineR)
library(ggbiplot)
##
## Attaching package: 'ggbiplot'
## The following object is masked from 'package:psych':
##
## reflect
birth.pca <- prcomp(corr_cl_data,
scale = T) # шкалирование нужно так как датафрейм подаем не шкалированный
summary(birth.pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 1.5936 1.0683 0.8825 0.8512 0.73273 0.52813
## Proportion of Variance 0.4233 0.1902 0.1298 0.1208 0.08948 0.04649
## Cumulative Proportion 0.4233 0.6135 0.7433 0.8640 0.95351 1.00000
# birth.pca$rotation
# графическая визуализация доли объясненной изменчивости для каждой из компонент
fviz_eig (birth.pca, addlabels = T, ylim = c (0,50))
В случае датасета birth алгоритм сработал достаточно хорошо, так как 3 компоненты объясняют 74% дисперcии (Cumulative Proportion of variance). 60% для первых 2 компонент означает что в исходных данных много скоррелированных колонок/ либо корреляция сильная.
contr_1 <- fviz_contrib(birth.pca, choice = "var", axes = 1)
contr_2 <- fviz_contrib(birth.pca, choice = "var", axes = 2) # axes - номер комоненты, из чего состоит каждая компонента?
contr_1
contr_2
pca_graph <- fviz_pca_var(birth.pca, col.var = "contrib")
pca_graph
На графике “Variables- PCA” выше визуализированы 1 и 2 главные
компоненты (по осям x и y, % показывает долю объясненной вариации
вариативности (дисперсии данных)). Длина стрелки (и цвет стрелки)
обозначает сколько корреляции было в переменной и каков вклад переменой
в дисперсию, объясненную PC1 и PC2. Видно три группы скоррелированы
переменных, причем hospstay и btw/gest отрицательно зависимы. Группа
bwt/gest образуют почти прямой угол, что говорит об их вероятной
нескоррелированности (по PC2).
На графике “Сontribution” видим из вклада каких переменных состоит компонента. Лучше, когда несколько переменных дают наибольший вклад в объясняемую компонентой дисперсию - в наших данных это наблюдается скорее для PC2 и PC3. Линия на графике показывает значения по оси Y, которое было бы в случае если все переменные вносили одинаковый вклад в компоненту.
library(ggbiplot)
biplot <- ggbiplot(birth.pca,
scale=0, alpha = 0.3, ellipse = TRUE, groups = cleaned_data$dead) +
labs(title = "PCA of yield contributing parameters", fill = "dead", color = "dead")
biplot
Высокие показатели для pltct/ apg1/ lowph возможно (направление для поиска) связаны с более благоприятным прогнозом (выживание, dead = 0).
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
ggplotly(biplot, textposition = "none")
ggplotly(biplot) %>% config(displayModeBar = F)
bipl <- ggbiplot(birth.pca, alpha = 0.2, scale=0, ellipse = TRUE, groups = cleaned_data$dead) +
geom_point(size = 2,alpha = 0.1, aes(col = cleaned_data$dead, text = paste0("ID: ", cleaned_data$id))) +
labs(title = "PCA of yield contributing parameters", fill = "dead", color = "dead")
## Warning in geom_point(size = 2, alpha = 0.1, aes(col = cleaned_data$dead, :
## Ignoring unknown aesthetics: text
labs(title = "PCA of yield contributing parameters", fill = "dead", color = "dead")
## $fill
## [1] "dead"
##
## $colour
## [1] "dead"
##
## $title
## [1] "PCA of yield contributing parameters"
##
## attr(,"class")
## [1] "labels"
ggplotly(bipl, tooltip = "text")
Прописано под графиками в п.7
Почему использовать колонку ‘dead’ для выводов об ассоциации с выживаемостью некорректно?
Так как 2 компоненты , что видно на графике, визуализируют только 2 компоненты, суммарно объясняя 61% дисперсии. Таким образом наблюдаемые тенденции могут служить лишь основанием для выбора дальнейшего направления поиска. Так как анализ эксплораторный, зависимости носят предположительный характер. При высоких значения dim можно бы ожидать что корреляция истинная.
UMAP позволяет проанализировать близость строк друг другу и плотнее сгруппировать разрозненные наблюдения на основании теории топологии, при снижении размерности. Из-за свойства сохранять локальные расстояния, метод подоходит для анализа отношния наблюдений/строк, но мало подходит для анализа соотношений колонок.
На графике UMAP просележиваются облака точек, однако при стандартных настройках облако исхода 1 (dead) полностью перекрывается с исходом 0, в то время как для PCA перекрытие облаков лишь частичное. PCA дает значительно больше информации, чем UMAP, за счет дополнительной информации о соотношении переменных
C использованием tidymodels:
library(tidymodels)
## ── Attaching packages ────────────────────────────────────── tidymodels 1.2.0 ──
## ✔ broom 1.0.6 ✔ rsample 1.2.1
## ✔ dials 1.3.0 ✔ tune 1.2.1
## ✔ infer 1.0.7 ✔ workflows 1.1.4
## ✔ modeldata 1.4.0 ✔ workflowsets 1.1.0
## ✔ parsnip 1.2.1 ✔ yardstick 1.3.1
## ✔ recipes 1.1.0
## ── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
## ✖ psych::%+%() masks ggplot2::%+%()
## ✖ scales::alpha() masks psych::alpha(), ggplot2::alpha()
## ✖ infer::chisq_test() masks rstatix::chisq_test()
## ✖ scales::discard() masks purrr::discard()
## ✖ plotly::filter() masks rstatix::filter(), dplyr::filter(), stats::filter()
## ✖ recipes::fixed() masks stringr::fixed()
## ✖ dials::get_n() masks rstatix::get_n()
## ✖ dplyr::lag() masks stats::lag()
## ✖ infer::prop_test() masks rstatix::prop_test()
## ✖ yardstick::spec() masks readr::spec()
## ✖ recipes::step() masks stats::step()
## ✖ infer::t_test() masks rstatix::t_test()
## • Learn how to get started at https://www.tidymodels.org/start/
library(embed)
umap_p <- recipe (~., data = corr_cl_data) %>% # "тех строка" подаем датасет
step_normalize(all_predictors()) %>% # нормализуем
step_umap(all_predictors()) %>% # проводим umap со станд настройками 15 и 0.01
prep() %>% # "тех строка"
juice() # приводим результаты umap к стандартизированному датасету
umap_p_0 <- umap_p %>%
ggplot(aes(UMAP1, UMAP2)) +
geom_point (aes(color = cleaned_data$dead),
alpha = 0.5, size = 2) +
labs(color = NULL) +
theme_minimal()
umap_p_0
Альтернативный способ, пакет (umap)
require (umap)
library (umap)
umap <- umap(corr_cl_data)
df <- data.frame(x = umap$layout[,1],
y = umap$layout[,2],
dead = cleaned_data$dead)
ggplot(df, aes(x, y, colour = dead)) +
geom_point()
make_umap <- function (corr_cl_data, x,y) {
um <- recipe (~., data = corr_cl_data) %>% # "тех строка" подаем датасет
step_normalize(all_predictors()) %>% # нормализуем
step_umap(all_predictors(), neighbors = x, min_dist = y) %>% # проводим umap
prep() %>% # "тех строка"
juice() # приводим результаты umap к стандартизированному датасету
um_p <- um %>%
ggplot(aes(UMAP1, UMAP2)) +
geom_point (aes(color = cleaned_data$dead),
alpha = 0.5, size = 2) +
labs(color = NULL) +
theme_minimal() +
geom_text(aes(x = 1, y = 3,
label = paste0("neigh: ",x,", min_dist: ", y)),
stat = "unique",
size = 6, color = "red")
um_p
}
make_umap (corr_cl_data, 3, 0.01)
#make_umap (corr_cl_data, 5, 0.1)
make_umap (corr_cl_data, 50, 0.001)
#make_umap (corr_cl_data, 50, 0.01)
С изменением параметров кол-ва соседей и минмального расстояния изменяется плотность точек, разреженность графика проекции. Самая “плотная” группировка достигается при уменьшении минимальной дистанции и колисечества соседей. Более равномерно распределенная проекуция при увеличении мин. дистанции и кол-ва соседей.
Сделаем новые колонки на основе данных анализированных ранее.
bwt_permut<- corr_cl_data #clone data for analysis
bwt_50 <- sample(corr_cl_data$bwt[1:246], 246, replace = FALSE) # сделали вектор на основе первой половину btw без повторов
bwt_permut$bwt_p50 <- c(bwt_50, corr_cl_data$bwt[247:492]) # # пермутировали 50% значений переменной btw без повторов, склеили
bwt_permut$bwt_p100 <- sample(corr_cl_data$bwt, 492, replace = FALSE) # пермутировали 100% значений переменной btw без повторов
PCA для данных с 50% и 100% пермутированными значениями bwt:
# PCA анализ
bwt_50 <- bwt_permut %>% select (!c(bwt, bwt_p100))
bwt_50.pca <- prcomp(bwt_50,
scale = T) # шкалирование нужно так как датафрейм подаем не шкалированный
summary(bwt_50.pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 1.486 1.0709 0.8876 0.8744 0.8090 0.66171
## Proportion of Variance 0.368 0.1911 0.1313 0.1274 0.1091 0.07298
## Cumulative Proportion 0.368 0.5592 0.6905 0.8179 0.9270 1.00000
biplot_50 <- ggbiplot(bwt_50.pca,
scale=0, alpha = 0.3, ellipse = TRUE, groups = cleaned_data$dead) +
labs(title = "PCA (50% permutated btw)", fill = "dead", color = "dead")
bwt_100 <- bwt_permut %>% select (!c(bwt, bwt_p50))
bwt_100.pca <- prcomp(bwt_100,
scale = T) # шкалирование нужно так как датафрейм подаем не шкалированный
summary(bwt_100.pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 1.3977 1.0451 1.0064 0.8695 0.8443 0.68746
## Proportion of Variance 0.3256 0.1820 0.1688 0.1260 0.1188 0.07877
## Cumulative Proportion 0.3256 0.5076 0.6764 0.8024 0.9212 1.00000
biplot_100 <- ggbiplot(bwt_100.pca,
scale=0, alpha = 0.3, ellipse = TRUE, groups = cleaned_data$dead) +
labs(title = "PCA (100% permutated btw)", fill = "dead", color = "dead")
biplot_50
biplot_100
ggarrange (biplot, biplot_50 ,biplot_100)
Кулумятивный процент объясненной вариации снижается (на примерно 5 и 10% для 50% и 100% пермутированных переменных соотвественно), что говорит об ухудшении качества модели. Наиболее значительно, ожидаемо, направление стрелки изменилось для самой переменной btw, особенно в случае 100% пермутации. Однако для других переменных, особенно с меньшей вероятностью скореллированных с btw, направление изменилось менее значительно (lowph) или практически не изменилось (hospstay).
Форма облаков и положение точек также претерпели некоторые изменения. По мере увеличения % пермутированнх значений- длина прочих стрелок меняется- они скорее всего теперь берут на себя больший вклад в долю общей обясненной вариации.
UMAP для данных с 50% и 100% пермутированными значениями bwt:
fun_umap <- function (data) {
umap <- recipe (~., data = data) %>% # "тех строка" подаем датасет
step_normalize(all_predictors()) %>% # нормализуем
step_umap(all_predictors()) %>% # проводим umap со станд настройками
prep() %>% # "тех строка"
juice() # приводим результаты umap к стандартизированному датасету
umap %>%
ggplot(aes(UMAP1, UMAP2)) +
geom_point (aes(color = cleaned_data$dead),
alpha = 0.5, size = 2) +
labs(color = NULL) +
theme_minimal()
}
fun_umap(bwt_50)
fun_umap(bwt_100)
#### 14. Импутация значений и иерархическая кластеризация
# датасет с импутацией значений вместо NA
birth_na_imp<- very_low_birthweight %>% select (hospstay, lowph, pltct, bwt, gest, apg1, dead) %>%
mutate_if(is.numeric, ~ (replace(., is.na(.), mean(., na.rm = TRUE)))) %>%
mutate (id = seq(1,671)) %>% relocate (id) %>%
select (!c(id, dead))
correlation_imp <- psych::corr.test(birth_na_imp,
method = "spearman", adjust = "fdr") # c поправкой на множ. сравнения
Визуализация взаимосвязей
# график матрицы корреляций
corrplot(corr = correlation_imp$r,
type="upper", order="hclust", method = "number", col = COL2('PiYG'),
p.mat = correlation_imp$p, sig.level = 0.01, insig = "blank")
# network (показывает силу связей и взаимоположение)
library(corrr)
net_imp <- birth_na_imp %>%
cor() %>%
network_plot(min_cor = .0)
ggarrange (net, net_imp)
Отмечается то, что пропала большая часть корреляци/взаимосвязи данных. в
частности достаточно сильная отрицательная взаимосвязь hospstay (времени
в клинике и переменных btw/gest)
Иерархическая кластеризацию на этом датафрейме.
birth_imp_scaled <- scale(birth_na_imp) # шкалирование
birth_imp_scaled_dist <- dist(birth_imp_scaled) # вычисление дистанций
df_dist_imp.hc <- hclust(d = birth_imp_scaled_dist,
method = "ward.D2")
fviz_dend(df_dist_imp.hc, horiz = TRUE,
cex = 0.6)
Одновременный график heatmap и иерархической кластеризации. Интерпретация результата.
library(pheatmap)
hierarchy_imp <- pheatmap(birth_imp_scaled,
show_rownames = FALSE,
clustering_distance_rows = birth_imp_scaled_dist,
clustering_method = "ward.D2",
cutree_rows = 5, #зададим количество кластеров для рядов
cutree_cols = length(colnames(birth_imp_scaled)),
angle_col = 45,
main = "Dendrograms for clustering rows and columns with heatmap")
hierarchy_imp
На графике видно как преобразование пропущенных зн методом замены на средние значения очень существенно исказило данные. При том что кластеры по переменным не значимых претерпели изменений, кластеры по рядам изменились существенно. Так же данные потеряли большую долю своей вариативности и усилилась тенденция к среднему, что выражается в большей “одноцветности” хитмепа по сравнению с графиком для датафрейма с удаленными значениями.
Как видно из графиков метод импутации хоть и позволят сохранить нам строки с частично пропущенными значениями, однако существенно искажает взаимосаязи между наблюдениями и может приводить к ложным направлениям поиска. Однако в части случаев, особенно когда данных мало и они невосполнимы, метод замены на среднее дает возможность сохранить информацию в случае неполных кейсов.
birth.imp.pca <- prcomp(birth_na_imp,
scale = T) # шкалирование нужно так как датафрейм подаем не шкалированный
summary(birth.imp.pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 1.5368 1.0158 0.9994 0.8468 0.8164 0.47350
## Proportion of Variance 0.3936 0.1720 0.1664 0.1195 0.1111 0.03737
## Cumulative Proportion 0.3936 0.5656 0.7320 0.8516 0.9626 1.00000
# birth.pca$rotation
# графическая визуализация доли объясненной изменчивости для каждой из компонент
fviz_eig (birth.imp.pca, addlabels = T, ylim = c (0,50))
pca_graph_imp <- fviz_pca_var(birth.imp.pca, col.var = "contrib")
pca_graph_imp
biplot_imp <- ggbiplot(birth.imp.pca,
scale=0, alpha = 0.3, ellipse = TRUE, groups = very_low_birthweight$dead) +
labs(title = "PCA (imputed NA)", fill = "dead", color = "dead")
biplot_imp
Процент объясненной дисперсии упал не значительно, для двух первых компонент с 61% до 57%, однако изменился паттерн взаимосвязей, который можно было бы предполагать исходя из графика, полученного в п. 7. Теперь hospstay вносит крайне малый вклад, и нет предполагаемого обратного характера взаимозависимости hospstay и gest/btw.
UMAP
umap_imp <- recipe (~., data = birth_na_imp) %>% # "тех строка" подаем датасет
step_normalize(all_predictors()) %>% # нормализуем
step_umap(all_predictors()) %>% # проводим umap со станд настройками 15 и 0.01
prep() %>% # "тех строка"
juice() # приводим результаты umap к стандартизированному датасету
umap_imp_p<- umap_imp %>%
ggplot(aes(UMAP1, UMAP2)) +
geom_point (aes(color = as.factor(very_low_birthweight$dead)),
alpha = 0.5, size = 2) +
labs(color = NULL) +
theme_minimal()
umap_imp_p
Большая часть данных с NA в переменных относилась к умершим детям,
поэтому при замене NA на средние точек, соотвествующих погибшим
(dead==1) стало значительно больше, и проекция UMAP визуально приобрела
большую “полярность” Групп 0 и 1.